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Abstract 

Information about the nature of the low-temperature anomalies and in particular the properties 
of the tunneling systems in silica and lithium silica glasses are revealed via computer simulations. 
The potential energy landscape of these systems is systematically explored for adjacent pairs of 
local minima which may act as double-well potentials (DWP) at low temperatures. Three different 
types of DWP are distinguished, related to perfectly coordinated silica, intrinsic silica defects, and 
extrinsic defects. Their properties like the spatial extension and the dipole moment are character- 
ized in detail. Furthermore, the absolute number of tunneling systems, i.e. symmetric DWP, is 
estimated. The results are compared with dielectric echo, specific heat and acoustic experiments 
on Suprasil I and Suprasil W. A semi-quantitative agreement for all relevant features is obtained. 
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INTRODUCTION 



For more than thirty years it is known that glasses display anomalous behavior for tem- 
peratures in the Kelvin regime. 1 For example, the temperature dependence of the specific 
heat turns out to be linear rather than cubic, as expected from the Debye model. The stan- 
dard tunneling model (STM) predicts many of these features. 2 ' 3 The key idea is to assume 
that the system can perform local transitions between adjacent configurations. Because of 
the low temperatures the crossing of the barriers occurs via tunneling. This feature can be 
characterized by a double-well potential (DWP) in configuration space, i.e. two adjacent 
local minima on the potential energy landscape (PEL). In general, the reaction coordinate 
corresponds to a highly cooperative process in which a group of neighboring atoms is in- 
volved. Among all DWP which are present in a given glassy configuration only those DWP 
with asymmetries less than, let's say, 2 K will contribute to the low-temperature anomalies. 
In contrast, typical energy scales for, e.g., Si0 2 are of the order of leV and thus many 
orders of magnitude higher. Thus one would expect that only a minor fraction of all DWP 
are relevant for the low-temperature anomalies. They are denoted two-level systems (TLS). 

There is no theory which predicts the properties of TLS from first principles without in- 
voking some model assumptions. There are, however, interesting approaches within different 
physical frameworks like a mean-field model 4 and random first order transition theory. 5 In 
the second approach a detailed picture of the PEL in terms of mosaics is derived. In their 
model the reaction coordinate of TLS involves the dynamics of O(10 2 ) molecules which only 
move a small fraction of a nearest-neighbor distance. Finally, in the STM some reasonable 
assumptions about the properties of the potential parameters of the TLS are made. 

Computer simulations are well suited to elucidate the properties of TLS. Two steps are 
involved. First, one has to generate a typical glassy configuration and, second, to identify 
nearby adjacent minima on the PEL and characterize the nature of the individual transitions. 
Qualitatively, one would expect that the active atoms for one transition are localized in 
some region in the glass. The first simulations used non-systematic search methods. 6,7 
In subsequent work Heuer and Silbey have introduced a method how to search TLS in a 
systematic way. 8 Applying this method to a binary mixture Lennard- Jones (BMLJ) system 
they were able to predict the real-space realization of the TLS 8 , the coupling to strain, 9 ' 10 
the prediction of tunneling properties in terms of material constants and, as a consequence, 
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a universal description of the low-temperature anomalies. 11,12 Important results about the 
properties of DWP in Lennard- Jones systems can be also found in Refs. 13_15 . For all 
simulations one has to take into account that the absolute number of TLS is so small that 
there is no chance to identify a sufficiently large number via simulations. In contrast, the 
number of DWP is very large. Therefore an important ingredient of these simulations is 
the formulation of a statistical method how to predict the properties of TLS from those 
of DWP. Using the significantly improved computer facilities during the last decade, the 
present authors have recently reanalyzed the BMLJ system . 16,17 In particular, it could be 
shown that (1) the systematic search algorithm indeed allows one to obtain an estimation of 
the absolute number of TLS, (2) the properties of the TLS basically do not depend on the 
energy of the configuration. Thus, the intrinsic computer limitations like the necessity of a 
very fast cooling procedure to generate the initial configuration or the choice of relatively 
small systems do not hamper the unbiased determination of TLS properties in the BMLJ 
system. 

A prototype glass-forming system is vitreous silica (Si0 2 ). Trachenko et al have charac- 
terized the TLS of silica in several respects and have shown that the microscopic nature is 
related to coupled rotations of SiC>4 tetrahedra. 18-21 A similar picture has been suggested 
from experimental neutron data 22 and from the analysis of localized vibrational modes. 23 
The transition events have been also analyzed via the activation-relaxation technique. 24 The 
present authors have conducted a systematic search for TLS in silica. 25 It has been possible 
to estimate the experimentally observed number of TLS which has turned out to be a factor 
of approx. 3 smaller than the experimental TLS density reported for Suprasil W. As argued 
in Ref. 25 there are two possible reasons for a slight underestimation by computer simulations 
which are connected with the estimation of the tunneling matrix element. Going beyond 
the WKB-approximation one might additionally take into account that (i) the saddles are 
broader (in terms of second derivatives) than the minima (see also Refs. 14,17 ) and (ii) there 
is a Franck-Condon factor, taking into account the relaxation of the phonon modes during 
the transition within a DWP. 17 Both effects would increase the estimation of the number of 
TLS and might at least partly account for the remaining difference. 

In general, one may distinguish intrinsic and extrinsic TLS. 26 The extrinsic TLS may be 
related to extrinsic defects like OH-impurities. Actually, comparing Suprasil I (1250 ppm 
OH-impurities) with Suprasil W (j 5 ppm OH-impurities) the experimentally determined 
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number of TLS is nearly twice as high. This may be related to the additional contribution of 
extrinsic TLS. 26 For the intrinsic TLS two different contributions may be of relevance. First, 
they may be related to TLS where intrinsic defects of the silica system are involved. They 
are characterized by non-tetrahedral local coordination like non-bridging oxygen atoms. 
Second, intrinsic TLS may also be present in defect-free silica configurations which recently 
have been analyzed in Ref. 25 . Intuitively, one might expect that defective configurations 
(intrinsic and extrinsic) are more efficient to locally reorganize and thus to form TLS. Both, 
the presence of intrinsic and extrinsic defects, will be relevant for a closer rationalization of 
the low-temperature anomalies in systems like Suprasil. 

In previous work on a BMLJ system the properties of extrinsic defects have been analyzed 
in detail. 27 For this purpose one has included a test Lennard- Jones particle with variable 
radius and analyzed the DWP, related to this particle. Some dramatic effects have been 
observed when analyzing a test particle with a radius which is e.g. 25% smaller than the 
radius of the smaller component of the BMLJ system (data compared with the properties 
of the majority component): (i) the probability for the formation of a DWP has increased 
by a factor of approx. 50, (ii) the number of particles, involved in the DWP, has decreased 
by more than a factor of 2, (iii) the average saddle height of the DWP has increased by a 
factor of nearly 2, (iv) the deformation potential has decreased by 25%. 

In this work we explicitly compare the nature of the TLS of defect-free silica with those of 
silica containing intrinsic defects and/or extrinsic defects, modelled by a small concentration 
of Li 2 0. It will turn out that all types of TLS have very different properties. On a qualitative 
level, many features for the extrinsic defects will turns out to be similar to what has been 
observed for the test particle in the BMLJ system, as sketched above. To underline the 
interpretation of the results, we compare them with the respective observations for the 
BMLJ system and with experimental data. 

TECHNICAL 

The molecular dynamics (MD) simulations have been conducted under NVE-conditions. 
The velocity Verlet scheme has been chosen to propagate the particles in time. The MD 
simulations were basically used to generate large numbers of independent structure, which 
then were analyzed for the occurrence of DWP. Periodic boundary conditions have been 
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applied. 



Binary mixture Lennard-Jones 

As a LJ-model glass former we chose a binary mixture system with 80% large A-particles 
and 20% small B-particles (BMLJ). 28-31 It is supposed to represent NiP (80% 62 Ni; 20% 31 P) 
but with a 20% higher particle density. 32 This system was first used by Kob and Anderson. 
The used potential is of the type 

V aP = 4 • e af3 [(a af3 /r) 12 - {a a p/rf\ + (a + b ■ r), (1) 

with <t ab = 0.8a AA , cr BB = 0.88a AA , e AB = 1.5e AA , e BB = 0.5e AA , m B = 0.5m A . A linear 
function a + b ■ r was added to ensure continuous energies and forces at the cutoff r c = 1.8. 
The units of length, mass and energy are a AA , m A , e AA , the time step within these units was 
set to 0.01. For the case of NiP the energy unit corresponds to 934 K and a AA is 2.2 A. The 
presented data are taken from our smallest simulated system with 65 particles. The finite 
size effects have already been analyzed in Ref. 16,17 and are small enough to be neglected. 

Silica 

The pure silica system has been modeled by the BKS-potential. 33 We have analyzed 
system sizes of N=150 and N=600 particles. If not mentioned otherwise the data in this 
work refer to = 150. The system has the standard density of 2.3 g/cm 3 and a short-range 
cutoff of the BKS-potential of 8.5 A. The starting configurations for the systematic search 
correspond to equilibrium configurations at 3000 K, which subsequently were minimized. We 
have obtained starting structures with and without intrinsic defects, i.e. deviations from a 
perfect tetrahedral coordination. Due to the fact, that defect related DWP may be relevant, 
though the number of defects is very small, 34 we do not only analyze the defect free DWP 
but also the defect related DWP. 

To define the coordination of a silicon atom in a classical pair potential it is necessary to 
introduce maximum bond cutoffs to determine if a bond exists or not. The minimum be- 
tween the first and the second nn-shell in the radial distribution functions for the minimized 
structures lies around 2.2A. This value has thus been used as a cutoff to define bonds, Si-0 
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distances below this value are considered bonds and larger distances are not. The general 
results do not critically depend on the used cutoff value as a broad minimum region between 
the two nearest-neighbor-shells exists. 

Lithium silicate 

For the alkali silica simulations a f 53 particle system with only 2 Li-atoms has been used. 
The low lithium concentration guarantees, that the Li-atoms behave almost independently 
and that the network structure is rather similar to that of pure silica. The potential has 
been taken from Habasaki and Okada. 35 Its alkali-free limit is close to the BKS-potential. 
The simulations have been conducted under the same conditions as for pure silica. 

Search for DWP 

The key idea of our search algorithm is to start from one minimized configuration, i.e. a 
local minimum of the PEL, and then to perform a specified number of MD steps. Afterwards 
the system is minimized again. When this procedure ends in a configuration which is different 
from the starting configuration it is checked whether there exists a saddle point between both 
configurations. If yes, one has identified one DWP. This procedure is repeated many times 
(typically: 100 times) in order to identify most if not all DWP. A DWP is characterized by 
three parameters: its asymmetry A, its potential height V and its distance between both 
minima d. More specifically we either use the Euclidean distance d or the mass-weighted 
distance d mwrp along the reaction path. 16 Since these differences are not relevant in the 
context of the present work they will be neglected in what follows. For the determination 
of the saddle a robust saddle search routine has been chosen. 31 

A DWP is kept for the further analysis if A < A max and d < d max where A max and d max 
are specified values. This limitation is essential to enable a systematic search. DWP with, 
e.g., very large values of the asymmetry are difficult to find such that a systematic search of 
all DWP (starting from a given configuration) is not possible within reasonable computer 
time. In contrast, if one is only interested in DWP within the range, specified above, one 
may hope to find (nearly) all DWP. As will be shown below, the TLS are confined to this 
range, i.e. TLS with d > d max are irrelevant. 
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A simple way to estimate the degree of completeness of the search procedure is to analyze 
how often a specific DWP is found during the repeated simulations. It has been shown that 
for appropriately chosen parameters a systematic search for DWP is indeed possible for the 
BMLJ system 16 as well as for silica 25 and properties of TLS, i.e. of nearly symmetric DWP, 
can be extracted. For a successful assessment of TLS properties one first has to deal with 
the strong positive correlations among all three parameters d, V, and A. Guided by the 
soft-potential model 36,37 we have mapped the distribution p(d, V, A) on a distribution of 
soft-potential parameters W2,w 3 ,W4 which fortunately to a very approximation turn out to 
be independent of each other. Based on the independent distribution functions Pi(wi) it is 
possible to generate a sufficiently large set of DWP with the same statistical properties as 
the initial set of DWP. In this way one can cover the full parameter range with high precision 
and in particular obtain the distribution of TLS as a subset of the whole parameter range. 
A closer description of this parametrization method as well as further technical details can 
be found in Refs. 16 ' 17 . 

The parameters, used for our simulations, are given in Tab.l. For the characterization of 
the DWP we first determine which particle moves most during the transition. This particle 
is denoted central particle. For BMLJ we have therefore distinguished A-type DWP where 
the central particle is an A-particle and, analogously, B-type DWP. For silica we have defined 
DWP ND as the subset of DWP where the initial and the final configuration are defect-free 
and no bond breaking occurs during the transition. The remaining set is denoted DWP WD . 
For lithium silicate we define DWP Ll as those DWP for which one of the two lithium atoms 
is the central particle. The number of DWP, found in our simulations, are also listed in 
Tab.l. 

RESULTS 

Properties of DWP parameters 

Before analyzing the total number of DWP, found in our simulations with the choice 
of A max and d max as listed in Tab HI one has to ask whether or not the search for DWP 
was complete. In Fig^we display the results for DWP ND and DWP WD when applying the 
completeness criterion. The curve for DWP ND behaves similarly to what we have found 
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# DWP 






BMLJ (A-type) 


521 


0.8 


0.5 


BMLJ (B-type) 


6001 


0.8 


0.5 


DWP ND 


250 


5 A 


1500 K 


DWP WD 


1419 


n 


r> 


DWP Li 


1885 







Table I: The number of DWP analyzed for all different systems (system sizes: N=65 for BMLJ; 
N=150 for silica; N=153 for lithium silica). They were recorded if their values of d and A satisfy 
d < d max and A < A max , respectively. 

for the BMLJ system. 16 It displays a maximum for a value much larger than one (here: 
approx. 30). Qualitatively, this means that whenever a DWP is present it is found quite 
often during the repeated runs. As a consequence it is statistically very unlikely that a DWP 
is not found at all. This is exactly the condition for a (nearly) complete search. This result 
has been already used in our previously published work. 25 In contrast, the probability that 
a DWP WD is found just once out of the 100 attempts is significantly larger than to be found 
twice or even more often. Thus, it is likely that several DWP WD are not found at all. A 
similar observation is made for the DWP Ll (data not shown). It is, however, unlikely that 
the true number of DWP is orders of magnitude larger than the number of DWP already 
found because otherwise the probability to find a DWP only once or twice would be much 
larger than 30% for DWP WD . Furthermore for the estimation of the TLS one mainly needs 
to have reliable information about DWP which are nearly symmetric and have a small value 
of d (see below for more details). It is likely that this subset of DWP is found with a higher 
probability because of its proximity on the PEL which further supports that at least the 
number of TLS, estimated below, is not far away from the true value. 

Based on the number of DWP, found in our simulations and listed in Tab [II and the 
number of analyzed initial configurations one can estimate the probability that a given 
configuration contains one DWP within the specified parameter range. For atomic systems 
like BMLJ it is useful to divide the number of DWP also by the number of particles per 
configuration to obtain a per-particle DWP probability. Thus, for the A- Type DWP we have 
defined the probability relative to the number of A-particles. The B-type DWP are treated 
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Figure 1: The histogram of how often a nearby minimum was found per starting minimum, using 
100 independent runs. 

analogously. For the lithium silicate system we consider the occurrence of DWP Ll relative 
to the two individual lithium atoms per configuration. For pure silica a SiC^-tetrahedron 
serves as an elementary unit because of its rigid character at low temperatures. Thus a silica 
configuration with N = 150 particle contains 50 individual elementary units. Actually, in 
some models of silica the tetrahedra are directly treated as rigid bodies. 18 We note in passing 
that the size of the elementary unit can be estimated from analyzing the material constants 
of silica. 38 Therefore, for a better comparison with other systems the number of DWP ND for 
defect-free configurations is related to the number of tetrahedra. It turns out that for a silica 
configuration with iV = 150 atoms the probability to find a DWP is 10 times larger if the 
configuration contains a defect. Thus the DWP can be to a good approximation exclusively 
related to the presence of a defect. On average, it turns out that for N = 150 a configuration 
contains approx. 1.2 independent defects. Using all these pieces of information one can 
express the number of DWP WD relative to the number of elementary units, which here are 
the number of independent defects. For all five different types of DWP the probability of 
DWP formation per elementary unit (within the specified range of parameters) is given in 
TablTTl Note that for DWP WD and DWP Ll we can only determine lower limits. 

For understanding the low-temperature anomalies one is interested in the number of TLS 
(here defined via A < 2 K) rather than in the number of DWP. In analogy to the DWP 
we define the TLS ND , TLS WD and TLS Ll . If the distribution of asymmetries were constant 
in the interval [0, A maa: ] the number of DWP with an asymmetry smaller 2K could be 
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DWP/elem.unit 
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TLS/elem.unit 


BMLJ (A-type) 


1.0 • 10" 3 


2.6 


1.1 • 10" 5 


BMLJ (B-type) 


4.0 • 10~ 2 


2.6 


5 • 10" 4 


DWP ND 


1.4 • 10~ 3 


2.0 


3.8 • 10" 6 


DWP WD 


> 5 • 10" 1 


1.1 


> 7.3 • 10~ 4 


DWP Li 


> 2.8 • 1Q- 1 


1.25 


> 4.7 • 10~ 4 



Table II: The probability that a specific elementary unit is the central unit for the transition of a 
DWP with asymmetry less than 2K, i.e. a TLS. /a characterizes the deviations from a constant 
distribution of asymmetries. 
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— DWP ND N=150 
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Figure 2: The distribution of asymmetries for the DWP ND and DWP WD . 

directly calculated by multiplying the above values with 2K/A ma:r . Closer inspection of the 
distribution shows, however, that the density is somewhat higher for very small asymmetries; 
see FigElfor DWP ND and DWP WD . The density for A < 300 K is roughly constant within 
statistical noise. One can define a factor /a which describes the increases of the density 
in the limit of low asymmetries as compared to the prediction for constant density in the 
whole interval. An estimation for this value is also shown in TablLTl As can be judged from 
FigO there is some statistical uncertainty of approx. 20%. The value of /a is particularly 
large for DWP ND and the DWP in the BMLJ system. Qualitatively, it is related to the 
barrier distribution in FigJHl Due to the strong correlation of asymmetry and barrier height 
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DWP very high barrier heights are very unlikely to have a small asymmetry. Thus a broad 
distribution of barrier heights implies that only a smaller fraction of DWP can be relevant 
for the low-temperature anomalies. Indeed, it turns out that to a good approximation the 
value of /a is proportional to the fraction of DWP with A < 1000 K in Fig 01 Now the 
number of TLS per elementary unit can be directly estimated. The results are also given in 
TabO 

As a main conclusion from this analysis it turns out that A-type DWP for the BMLJ 
system and DWP ND are very rare. The other extreme are the extrinsic lithium defects 
DWP Ll as well as the intrinsic defects DWP WD in silica. Thus the presence of a small 
particle in a disordered network or the breaking of this network, e.g. via non-bridging 
oxygens (see below), prompts the formation of bistable modes, i.e. DWP. Hence via local 
motion of the defect the system can acquire a new metastable position. Also the probability 
of the formation of B-type DWP for BMLJ is dramatically enhanced. This is consistent 
with the view that the B-particles can be regarded as defects among the majority of the 
A-particles. We note in passing that for the lithium silicate system also the number of DWP 
with silicon or oxygen as central particles is increased by approx. 50% if compared to the 
number of DWP for an average silica configuration. This may be related to the fact that 
the presence of lithium automatically implies a breaking of the silicate network and then 
the additional intrinsic defects may give rise to an increased number of DWP. 

A central DWP parameter is the barrier height. In Fig 01 we display the distribution 
of barriers for the different types of DWP. Interestingly, the absolute values of the barrier 
heights are very different. For DWP ND the transition between both minima can be described 
as a coupled rotation of several Si0 4 -tetrahedra. 25 Obviously, this can be achieved by sur- 
mounting only a relatively small energy barrier. A rationalization for the higher barriers for 
dwp wd and DW P Li will be presented further below. In Ref. 39 the distribution for DWP ND 
has been estimated to be confined to approx. V < 1000 K. This is in very good agreement 
with the data shown in FigO 

Furthermore, one may ask whether the properties of the DWP depend on the energy of 
the starting configuration. Therefore we analyze the average barrier height of the DWP in 
dependence of this initial energy; see Fig0J The energies correspond to the relevant energies 
at low temperatures. Actually, the previous analysis about the relaxation properties of 
silica has revealed that the PEL of silica has a low-energy cutoff which for the present 



11 




V/K 

Figure 3: Barrier distribution for the different types of DWP. 

system size can be estimated to be around -2895 eV. 40 For all three systems there is no 
significant dependence of the average potential height on energy relative to the level of 
statistical uncertainties. Thus, one may conclude in analogy to our previous results for the 
BMLJ system, that the fast cooling rate in computer simulations is no serious problem for 
the quantitative analysis of DWP. Also different quantities, analyzed along this line, do not 
show a significant energy-dependence. 

Finally, we have analyzed the distance d between the two minima of the DWP. The 
results are shown in Fig El Nearly all DWP ND display distances between 1.5 A and 4 A. 
This implies that the motional pattern of coupled rotations of tetrahedra requires, on the one 
hand, some minimum length and, on the other hand, does not extend beyond some upper 
limit. In contrast, in particular for the defect dynamics, i.e. DWP WD a larger variance 
of possible distances is observed. This hints towards a larger number of different motional 
mechanisms in the presence of defects. 

Microscopic properties of DWP 

A central question deals with the number of atoms involved in the transition. This can 
be captured by the participation ratio. In literature different definitions can be found. Here 
we use R H and R as defined in Tab Mil d max is the distance moved by the central particle 
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Figure 4: Dependence of the average DWP barrier on the energy of the initial configuration. The 
DWP Li data are shifted by -200 eV to handle the different minimum energies. 





R H 


R s 


BMLJ {A-type) 


6.7 


16.1 


BMLJ (B-type) 


3.0 


9.0 


DWP ND , N = 150 


9.1 


24 


DWP ND , N = 600 


9.0 


31 


DWP WD , N = 150 


6.4 


19 


DWP Li , N = 153 


1.9 


3.6 



Table III: Participation ratios for the different systems studied in this work. 

and di the distance moved by the i-th particle. The results are listed in Tab II I II 

As the most prominent feature the transition in DWP Ll is extremely localized. Thus 
one may speak of a single-particle transition. This agrees with our previous observation for 
the BMLJ system. There the B-type DWP, corresponding to the dynamics of the small 
B-particles, are localized, too. In contrast, larger participation ratios can be found for the 
A-type DWP in the BMLJ system as well as the DWP in the silica system. Similar values 
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Figure 5: The distribution of distances for different DWP. 

have been obtained in Refs. 19 ' 41 . Taking into account that for silica the elementary units are 
SiC>4 tetrahedra and not individual atoms, these values should be considered as upper limits 
for the number of independent degrees of freedom, involved in the transitions. For DWP ND 
we have performed a systematic search for N = 150 and N = 600 particles. In agreement 
with other observables finite size effects are only weak. Thus, one may conclude that neither 
the finite size of our sample nor the fast quenching rate, used in computer simulations, have 
a relevant effect on the properties of the DWP. 

For the case of the BMLJ system the participation ratio is slightly lower for DWP with 
smaller distances d (unpublished work). Since TLS have smaller distances than the average 
DWP, the above values for the participation ratio should be regarded as upper limits for 
TLS. Thus, one may conclude that in all cases one has less than, let's say, 10 particles which 
participate in a tunneling mode. 

For a closer analysis of the nature of the transitions we have analyzed whether the transi- 
tion between both minima for the different particles is basically a straight line in real space 
or whether it is strongly curved. For this study we compare for all particles the transition 
vector from the first minimum to the saddle with the corresponding vector from the sad- 
dle to the second minimum and determine the angle a between both vectors. In case of 
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Figure 6: The average cosine of the angle between the transitions from the first minimum to the 
saddle and the transition from saddle to the second minimum. The particles are sorted according to 
their absolute value of the Euclidean distance when moving between both minimum configurations. 

a straight line both vectors are parallel, i.e. a = 0. In the other extreme case a particle 
performs a forward-backward motion such that both vectors are antiparallel, i.e. a = tt. 
Furthermore, we sort the particles according to the Euclidean distance between the initial 
and final position. The particle with index i — thus corresponds to the central particle, 
the particle with index i — 1 to the particles moving second most and so on. This allows us 
to average over all DWP. The results are shown in FigEl 

Clearly, for all systems the central particle of a DWP moves along a nearly straight line 
during the transition ((cos a) ~ 0.9). However, for the other particles major differences exist. 
For A-type DWP of the BMLJ system all particles perform a relatively straight transition 
path. This reflects the cooperative nature of the transition. In contrast, for B-type DWP 
already the transition path for the particle with i = 1 displays a strong curvature. This 
is compatible with the single-particle nature of the DWP. During the transition of a small 
B-particle between two local minima of the PEL the surrounding A-particles mainly retreat 
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during the transition thereby reducing the value of the saddle energy and finally go back to 
a similar position as before. 

Not surprisingly the same difference is observed when comparing DWP ND and DWP WD 
with DWP Ll . The one-particle type transitions for DWP Ll are, in analogy to the B-type 
DWP, connected with a strong curvature for the transition of the remaining particles. There 
are, however, additional differences between DWP ND and DWP WD for small particle index i. 
The coupled tetrahedra-rotation for DWP ND implies a similar behavior of the most mobile 
particles. In contrast, for the transitions involving defects the degree of cooperativity is 
somewhat smaller. 

Comparing DWP ND with the A-type DWP from the BMLJ system it turns out that that 
the angles are generally larger in the silica system. The motion is thus more curved, which 
seems reasonable due to the network structure and fits to displacements due to rotating 
tetrahedra. 

Microscopic mechanism of the DWP transition 

In Ref. 25 it has been shown that for DWP ND the crucial step during the transition is a 
Si-O-Si bond flip. Correspondingly, in more than 99% of all cases the central particle is 
an oxygen and only starting from particle index % — 10 (see above for its definition) the 
probability to find a silicon atom is close to the statistical value of 1/3. Qualitatively, the 
dominance of oxygen atoms as the most displaced particle is compatible with the picture 
of cooperative tetrahedra rotations which mainly involve the dynamics of oxygen atoms. 
Interestingly, it turns out that oxygen atoms, acting as central atoms, possess a significantly 
shorter Si-0 bond length. This structural motif seems to enhance the probability for the 
formation of a DWP. 25 

No simple picture for DWP WD can be formulated. Rather we have observed a variety of 
different mechanisms from a close inspection of 30 DWP. Most of the times the transition 
could be described by one of the following three mechanisms: (1) A bond between a silicon 
and a threefold coordinated oxygen breaks, the silicon forms a new bond with another 
oxygen. This oxygen is now threefold coordinated. (2) A bond with a threefold coordinated 
oxygen breaks and a silicon center changes its coordination from four to three or from five 
to four. (3) A dangling oxygen (coordinated to only one silicon) forms a bond to a three- or 
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Figure 7: Partial pair correlation function <?Li-o( r ) f° r the central lithium of the DWP Ll in 
comparison with the average radial distribution function. Only the first peak is shown as no 
differences exist in other sections. 

fourfold coordinated silicon resulting in a four- or fivefold coordinated silicon and sometimes 
another oxygen from that silicon becomes either a dangling oxygen itself or coordinates to 
another silicon. The variety of these scenarios in part also reflects the different types of 
DWP. Not surprisingly, in this case the probability that the central particle is an oxygen 
atom is only slightly enhanced as compared to the statistical value (75%). 

To study the case of DWP Ll we have analyzed the partial Li-0 pair correlation function 
for the minimum configuration. The results for the nearest-neighbor shell are shown in 
Fig|7| The first peak of this double-peak structure is mainly related to non-bridging oxygen 
the second peak to bridging oxygen. The density of oxygen atoms, related to the first peak 
(and thus being mainly non-bridging oxygen), is decreased by approx. 25% when a DWP 
is present. The integration of both peaks shows that density is transferred from the first 
to the second peak, but no density is lost from the first two peaks. Thus a reduction of 
the density of non-bridging oxygen and a corresponding increase of the number of bridging 
oxygen is a prerequisite for the formation of a DWP. Interestingly, the average number of 
oxygen atoms in the nearest-neighbor shell decreases from 3.9 to 3.4 when comparing the 
minimum with the saddle. Thus, during the transition the lithium ion changes part of its 
oxygen neighborhood. 
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Dipole Moment 

Of major experimental interest is the coupling of the TLS to electric fields via their dipole 
moment. It can be easily determined via 

M = J2<lidi (2) 

% 

—* 

where denotes the partial charge and di the translational vector of particle i during the 
transition. Of interest is the absolute value of the dipole moment M = J \M\ 2 which can be 
calculated in a straightforward manner from Eq|21 More precisely, we have averag ed M 2 and 
finally calculated the square root. The partial charges in our potentials are 0.87e for lithium, 
2.4e for silicon and -1.2e for oxygen. Note that a single-particle transition with d = lA and 
unit charge would correspond to 4.8 Debye. The results are shown in Fig|K]for the different 
types of DWP. We have averaged DWP with similar distances d, in order to keep track of a 
possible dependence on d. To a good approximation the data can be described by a linear 
relation, i.e. 

M = (qd (3) 

where q is the average partial charge and ( a dimensionless proportionality constant. In 
order to convey a feeling for the absolute value of M it may be instructive to estimate M for 
the simple scenario for which all particles move independently, i.e. the di are uncorrelated. 
In this limit one obtains the same expression as Eq|3]with £ = 1. Thus the proportionality 
to d is a generic property which results from the definition of the dipole moment whereas 
the value of ( contains information about the motional mechanism. 

For DWP Ll one would expect M(2A) m 8D (using the partial charge of lithium) which 
is close to the numerically found value. Thus we find £ ~ 1. In any event, for a localized 
transition any correlation effects between adjacent particles cannot be relevant for the dipole 
moment. The other extreme case is DWP ND for which the dipole moment is approx. 10 
times smaller than expected for the uncorrelated case. This clearly shows that the dynamics 
is highly cooperative. In general, two effects may give rise to the reduction. First, particles 
of opposite charge move in the same direction; second, identical particles move in opposite 
directions. The latter case is relevant for DWP ND because of the dominance of tetrahedra 
rotations, involving significant oxygen motion. In case of ideal rotations around a threefold 
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Figure 8: The dipole moment for DWP ND DWP WD and DWP Ll in dependence of the distance 
between both minima. Included is a line of slope 1 and the average distances, estimated for the 
TLS. 

axis the resulting dipole moment would be zero. Due to the non-symmetric nature of the 
actual rotational axis some residual dipole moment remains. 25 

For DWP WD the reduction as compared to the statistical value is only a factor of approx. 
3 (taking into account the somewhat higher effective charge because also silicon atoms 
can move significantly). Despite the similarity in the participation ratios of DWP ND and 
DWP WD this shows that the nature of the cooperativity is significantly different. Beyond 
the rotation of tetrahedra it is the transfer of the defect structure which plays an essential 
role and which does not possess the cancelation effect of pure tetrahedral rotations. 

For comparison with experimentally determined dipole moments one has to take into 
account the dependence on the distance between both minima of the DWP. As already 
discussed above the distance for TLS is much smaller than that for the average DWP found in 
our simulations. This is a direct consequence of the generic correlations between asymmetry 
and distance. Using the parametrization method, sketched above, it is possible to calculate 
the distribution of TLS from the distribution of DWP. For the present analysis we are more 
specific with the definition of TLS. Beyond a near-degeneracy of the lowest two eigenstates 
(energy difference < 2 K) of the TLS we also require that the relaxation time r is smaller 
than 5s. 17 This means that the TLS is active on the time scale of typical experiments. 
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Figure 9: The distance distribution for a subset of DWP fulfilling different criteria as mentioned 
in the inset. 





d-TLS 


C 


M 


TLS ND 


1.20 A 


0.1 


0.65 D 


xls wd 


0.93 A 


0.3 


3.5 D 


TLS Li 


0.83 A 


1.0 


5.0 D 



Table IV: The dipole moments for the different types of TLS. Included are the reduction factors 
£ due to correlated motion of the participating atoms. 

This additional criterion in particular excludes DWP with very high barriers and/or large 
distances. This reduces the values of d for the relevant TLS even further. The results for 
D WP ND are shown in FigEl 

Using the parametrization method we first have recalculated the distribution of DWP, 
compatible with our limiting values d max , A max . This is the right curve in Fig El It is, 
of course, very similar to the original density distribution in FigEl Using the additional 
restriction E < 2 K one obtains a set of DWP which is shifted towards smaller values of d. 
The final set of DWP with the additional restriction r < 5s basically excludes DWP with 
d > 2 A. The average distance of these TLS is drLS = 1.2A. Thus together with FigJElone 
may estimate a dipole moment of 0.65 D (see Tab llV)) for TLS in the non-defect case which 
is close to the value of 0.6 D, reported in Ref. 26 . 
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We have performed the same analysis for the other two cases. The results for the values 
of dns and the dipole moment are given in Tab II VI Interestingly, the experimental value 
for OH-defects is approx. 4 D which is close to the value for TLS Ll as well as TLS WD . Thus, 
it appears that OH-defects and Li-defects have a similar transition mechanism during their 
tunneling motion. 

In recent work Eq|3]has been derived, using a somewhat different notation. 42 The deriva- 
tion was based on the ideas of the random first order transition theory. The parameter ( 
was empirically introduced as a phenomenological constant (( ~ 0.1) which is needed to 
recover the experimentally determined dipole moment of TLS ND . It has been postulated 
that (q is the effective partial charge, stemming from a Coulomb charge on a bead, which 
is much smaller than the electron charge e. 42 In contrast, in the present work we can show 
that the smallness of ( can be fully explained by the complex motional mechanism of the 
TLS transitions in silica. 



DISCUSSION AND SUMMARY 



Via extended computer simulations we have analyzed the properties of the DWP in silica 
and lithium silicate. Most importantly, the typical limitations of computer simulations, 
i.e. finite-time (resulting in fast quenching rates) and finite-size effects, do not hamper 
the characterization of the DWP within statistical uncertainties. The latter point can be 
concluded from the fact that the DWP properties do not depend on the energy of the 
initial glassy configuration and, in general, for small systems the energy is most relevant 
for predicting the dynamic properties. 31 This is analogous to the case of BMLJ. We should 
note in passing that for the BMLJ system properties like the vibrational density depend 
on the energy of the configuration and thus on the cooling rate. 43 In contrast, for silica 
even the vibrational properties are to a large extent independent of energy. 44 A quantitative 
estimation of the number of DWP was possible for the BMLJ system and the DWP ND . 
In contrast, for DWP WD and DWP Ll a complete search for DWP was not possible within 
accessible computer times. However, one could obtain reasonable lower bounds for the 
number of DWP. 

The lithium ions in lithium silicate may be regarded as a prototype system for extrinsic 
defects, immersed in a pure glass former. The DWP with lithium as central particle, i.e. the 
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TLS mainly correspond to single lithium transitions. On a qualitative level the situation 
is similar to the case of the small B-particles in the BMLJ system which may be regarded as 
defects being added to the larger A-particles. Also in this case the average DWP resembles a 
single-particle transition (this holds even more for a small test particle in the BMLJ system). 
Furthermore, in both cases the probability to form DWP is significantly higher than in the 
remaining glass (A-particles in BMLJ and silica in lithium silicate). Thus extrinsic defects 
are very efficient in prompting the formation of DWP. 

The transitions in DWP ND and DWP WD as well as the A-type DWP in the BMLJ system 
involve the displacement of a somewhat larger number of particles as reflected by larger 
values of the participation ratio. This spatial extension, however, is still relatively small as 
compared to some general predictions about the nature of TLS. 5 

A closer analysis shows that DWP ND and DWP WD , i.e. DWP in a configuration without 
defects and with defects, respectively, behave very differently. The transition dynamics for 
DWP ND can be characterized as coupled rotations of Si0 4 -tetrahedra. In particular this 
gives rise to a very small value of the dipole moment of this transition, i.e. a very small 
effective charge transport. In contrast, for DWP WD the dynamics is largely determined by 
defect-specific modes which effectively give rise to a much larger dipole moment. One may 
speculate that this is the reason why on average the potential barrier for DWP WD is much 
larger than that of DWP ND . A significant charge transport involves a major reorganization 
of the whole system which typically will be connected with a significant activation energy 
due to the dominant character of the Coulomb energy. 

To which degree do our simulations explain the results of low-temperature experiments 
for vitreous silica? 

Defect-free TLS: As already mentioned above we find approx. 4 • 10~ 6 TLS /tetrahedron 
without defects. As discussed in Ref. 25 this translates into an effective density of TLS P, ac- 
cessible from acoustic experiments, which within a factor of 3 agrees with the experimentally 
observed value for Suprasil W. 26 

Dipole moments: The dipole moment for TLS ND is close to the experimental value, ob- 
tained for Suprasil W. 26 Increasing the concentration of OH, i.e. using Suprasil I, one 
obtains a second contribution with a dipole moment which is larger by nearly one order of 
magnitude. 26 Naturally, it has been related to the tunneling of the OH-impurities. Inter- 
estingly, similar dipole moments are observed for TLS Ll and TLS WD . Since TLS with such 
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large dipole moments are absent in Suprasil W one may conclude that TLS WD do not play 
a major role. 

Relevance of intrinsic defects: The previous conclusion can be directly checked by ana- 
lyzing the estimated number of TLS WD . It would be of the same size as the TLS ND from the 
undistorted silica network if there exists at least one intrinsic defect per 200 SiO^tetrahedra 
(please have in mind the uncertainty in the determination of the number of DWP WD and 
DWP Ll as discussed above). Analysis of molecular dynamics simulations of silica, using the 
metabasin approach, suggests that around T g there is roughly 1 silicon defect per 300 tetra- 
hedra (data not shown). A different way of extrapolation for BKS-data yields a significantly 
smaller fraction of defects. 34 Furthermore, density functional calculations in general lead to 
a smaller number of defects than using the BKS-potential. 45 To the best of our knowledge 
no definite experimental information is available for the number of defects in pure silica. 
Thus, combining our results for the probability of generating TLS WD with the (somewhat 
uncertain) absolute number of intrinsic defects, the numerical results are at least consistent 
with the experimentally observed absence of TLS WD (using the conclusion from above). 

Relevance of extrinsic defects: The number of TLS Ll would be close to that of TLS ND if 
there is one lithium ion per 125 tetrahedra. For Suprasil I (1200 ppm OH-defects by weight) 
one thus has one OH-molecule per 240 SiCVtetrahedra. Experimentally it has been found 
via dielectric echo experiments 26 as well as specific heat experiments 46 that the number of 
extrinsic TLS is half the number of intrinsic TLS. Thus an equal contribution of extrinsic 
and intrinsic TLS would require one OH-molecule per 240/2=120 tetrahedra which is close 
to the value we obtained for an equal contribution of TLS Ll to TLS ND . Of course, this 
perfect agreement is to some extent accidental, given the different nature of lithium and 
OH-defects and the non-completeness of the search. In any event, the order of magnitude 
seems to be fully compatible. 

Coupling to acoustic modes: It is estimated from experiments that the deformation po- 
tential is approx. 40% smaller for extrinsic TLS than for intrinsic TLS (more specifically, 
TLS ND , according to the discussion above). This result is obtained, first, from the relaxation 
of electric echoes 26 and, second, from the comparison of the specific heat and the thermal 
conductivity for Suprasil I and Suprasil W. The difference between both probes is much 
smaller for the thermal conductivity. 47 It is known from theoretical considerations that the 
deformation potential scales with the average distance d. 9 The different values of d in Tab II VI 
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for TLS as compared to TLS suggest that indeed the deformation potential for extrinsic 
TLS may be significantly smaller than for intrinsic TLS. 

In summary, present-day computer simulations are able to reveal many microscopic prop- 
erties of two-level systems in glasses in the Kelvin regime and enable a quantitative com- 
parison with experimental data. Via the complementary information from theory and ex- 
periments a detailed knowledge about the underlying nature of tunneling systems becomes 
accessible. In particular the magnetic field dependence of polarization echo experiments 
may be promising to yield further insight about the microscopic nature of TLS from the 
experimental side . 48 > 49 In any case, based on the additional information from this type of 
simulations the tunneling systems need no longer be considered as phenomenological entities. 
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